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Abstract 

We incorporate hydrodynamic interactions in a structure-based model of ubiquitin and demon- 
strate that the hydrodynamic coupling may reduce the peak force when stretching the protein at 
constant speed, especially at larger speeds. Hydrodynamic interactions are also shown to facilitate 
unfolding at constant force and inhibit stretching by fluid flows. 



I. INTRODUCTION 



It has been widely recognized that the water environment affects the energy landscape 
and functionality of biomolecules in a profound way. There is, however, another solvent- 
related effect that is considered less frequently: the hydrodynamic interactions (HI) between 
individual segments of a biomolecule that moves. These interactions may affect the dynamics 
of conformational changes because any motion of one segment generates a local fluid flow 
which influences another segment. 

The presence of HI is known to affects dynamic properties of soft mater. For instance, HI 
modify the values of diffusion coefficients in colloidal suspensions^, affect the characteristics 
of the coil-stretch transition in polymers^, change the kinetic pathways of phase separation 
in binary mixtures^, and alter the kinetics of macromolecule adsorption on surfaces^. Much 
less is known about the role of HI in protein folding and unfolding processes. DickinsorP and 
TanakeP speculated that HI might affect the kinetics of protein folding, but the actual nu- 
merical assessment of the role of HI has come with the paper by Baumketner and HiwatarP. 
They have considered coarse grained models and found that HI delay folding of a /3 hairpin 
but do not affect folding of the a-helix. 

In this paper, we consider mechanical stretching of proteins and study the relevance of 
HI to the process. The stretching can be accomplished in several ways and we discuss three 
modes: at constant speed, at constant force, and through fluid flow. We chose ubiquitin as a 
model system, since there is a large body of experimenta l 8 * 9 * 1 ^ *^ and theoretica l 12 * 13 * 14 * 15 * 16 * 17 * 
data on its unfolding. A coarse-grained, Go-type modeP^ of a protein is used, constructed 
based on the knowledge of its native state. The Go models have been shown to give sur- 
prisingly good agreement with both the experimental results^* 19 * 20 * and all-atom molecular 
dynamic simulations^ when it comes to stretching. 

We outline the model in section 2, then introduce two different ways of tracking the 
evolution of the system: through Langevin Dynamics and Brownian Dynamics in sections 3 
and 4 respectively. In the following sections we discuss results pertaining to the three modes 
of stretching and demonstrate that HI can take many roles: they inhibit unfolding by fluid 
flow, but make the constant force stretching faster. At constant speed, they reduce the peak 
force if the speed is sufficiently high. This Hl-related reduction in force may be downplayed 
in the all-atom simulations of titin by Lu and SchulterP^ which would provide part of an 
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explanation for the excessively large forces obtained in these studies. 



II. THE COARSE-GRAINED PROTEIN MODEL 



In our simulations, we use the coarse-grained, Go-type model of a protein. In the model, 
each residue is represented by a single bead centered on the position of the C a atom. The 
successive beads along the backbone are tethered by harmonic potentials with a minimum at 
3.8 A. The other interactions between the residues are split into two classes: native and non- 
native. This determination is made by checking for native overlaps of all atoms in aminoacids 
when represented by enlarged van der Waals spheres as proposed in reference^. The amino 
acids, i and j, that do overlap in this sense are endowed with the effective Lennard- Jones 
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potential Vij = 4e ~ [y 1 ) ■ The l en gth parameters <7jj are chosen so that the 

potential minima correspond, pair-by-pair, to the experimentally established native distances 
between the respective aminoacids. In order to prevent emergence of entanglements, the 
non-native contacts are endowed with a hard core repulsion described by the r~j 2 part of 
the Lennard- Jones potential combined with a constant shift term that makes the potential 
vanish smoothly at 4 A. The specificity of a protein is contained in the length parameters cr^. 
The energy parameter, e, is taken to be uniform. We take e/ks = 900A", which correlates 
well with the data on titin and ubiquitin unfoldin g ^ 14 * 23 !. Thus the reduced temperature, 
T = kgT/e of 0.3 should be close to room temperature. All of the simulations reported 
here were performed at this temperature. Various simulation methods to study the dynamics 
of the system are outlined in the following sections. 



III. LANGEVIN DYNAMICS (LD) 

In this case, the dynamics of a protein is assumed to be governed by the Langevin equation 

mr\ = - 7 (r< - u(r<)) + F c { + T { . (1) 

Here, ri is the position of the z'th aminoacid, is the net force on it due to contact 
potentials, 7 is the friction coefficient, and 11(17) denotes the solvent flow field. Finally, T is 
a white noise term with the dispersion obeying 

(T i (t)T j (t')) = 2k B T' r 8(t-t')I5 ij 
3 



where I is the identity matrix. The white noise term mimicks the effect of the random 
collisions of the aminoacids with the surrounding solvent at the same time serving as a 
thermostat of the system. However, this scheme completely neglect the effects of HI which 
may exist in a real system, when the motion of one particle induces the flow influencing the 
dynamics of all the other particles. 

In the simulations, the friction coefficient 7 is taken to be equal to 2m /r where r = 
\jma 2 /e ~ 3ps is the characteristic time scale of oscillations in the Lennard- Jones well. The 
parameter cr=5 A used in the above definition is a characteristic value of in the system. 
The selected value of 7 corresponds to a situation in which the inertial effects are smalP^^ 
but the damping action is not yet as strong as in water. The equations of motion are solved 
by a fifth order predictor-corrector scheme. 

Let us note, however, that although Langevin dynamics is commonly used in simulations 
of biological systems, the validity of this approach is not always well-established. Namely, 
as already noted by Lorentz^, the Langevin equation in the above form may only be used 
if there is a separation of time scales between the relaxation time of the particle (i.e. the 
bead) velocity t v = — = and the viscous relaxation time of the solvent, t v = a 2 ^, 
where p is the density of the particle, a its radius, and p s - the density of the solvent (see 
also the thorough discussion in the book by MazcP^.) The ratio of those two time scales is 
proportional to p/p s . Since the densities of proteins are only about 50% higher than those 
of the surrounding liquid-^, there is no separation of time scales between the relaxation of 
fluid variables and those of the bead and, strictly speaking, instead of Eq. Q one should use 
the generalized Langevin Equation involving a memory kernel 

mT\(t) = - /"dt , e(t-0(ri(0-u(r i ,0)+F. ? (*)+r i (t) • (2) 
Jo 

where the noise is again Gaussian and related to the dissipative term through the generalized 
fluctuation-dissipation relation 

(T i (t)T j (t')) = k B TC(t-t')i8 ij . 

Such an approach is naturally much harder to implement (see, however,^) and thus the 
ordinary Langevin description as in Eq.Q is usually resorted to. In this paper, we show 
that in protein unfolding simulations the results of a simple Langevin Dynamics ([I]) are 
consistent with those of Brownian Dynamics (see Sec. 4). The latter is not affected by the 
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solvent and particle inertia effects, hence the agreement between the two methods seems to 
imply that non- instantaneous response of the solvent to the change of particle velocity does 
not play any important role in the protein unfolding processes. 

IV. BROWNIAN DYNAMICS (BD) 

If the momentum relaxation time scale (r„) is small in comparison to the time scales char- 
acterizing the conformational evolution of the system (r c ), it is appropriate to describe the 
dynamics in terms of equilibration of the particle configurations only. The exact definition 
of r c is problem dependent: if the protein is stretched by the flow U then r c = fj, if a force 
F acts on the molecule then r c = and if the Brownian diffusion plays the central role in 
particle evolution then r c = a 2 /D, (where D is the diffusion constant and a the radius of 
the bead). The algorithm for simulations of the evolution of particle positions in this time 
regime has been devised by Ermak and McCammoiP. The displacement of particle % in 
time step At (in the absence of the flow) obeys 

n -*°i= E(Vi " Dj)At + ^ E D° • FjAt + B h (3) 

where the index denotes the values of respective quantities at the beginning of the time 
step, Fj is the force exerted on particle j by other particles, D is the diffusion tensor and 
B - a random displacement given by a Gaussian distribution with an average value of zero 
and covariance obeying 

< BjBj >= 2D?-At (4) 

It is nontrivial to generalize the above expression to incorporate the effects of a general 
external flow fielcP^. However, in the present case, we will be interested only in the uniform 
flow, in which case one gets simply 

r, - r° = UAt + E(V, • D° ) At + ^ E ■ F°At + B„ (5) 

where U is the flow velocity. If the diffusion tensor is nondiagonal, there exists a coupling 
between the force acting on the particle j and the displacement of particle i (cf. Eqj3|. 
This coupling, mediated by the solvent, is commonly referred to as the "hydrodynamic 
interactions" . 
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Note that without the HI, the diffusion tensor is simply 

k B T. 
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(6) 



and we recover the overdamped limit of Eq.(l). The diffusion tensor D depends in a com- 
plicated nonlinear way on the instantaneous positions of all particles in the system. For a 
system of spheres, exact explicit expressions for the diffusion tensor exist in the form of 
the power series in interparticle distances, which may be incorporated into the simulation 
scheme^- 1 | 32 | 33 | 34 | 35 | 36 | S7 1 . Here we adopt a pairwise, far-field approximation of D proposed by 
Rotne, Prager^ and YamakawaP*' 
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where r, 



rj and a represents the hydrodynamic radius of a bead. Since the above 



expression is exact only in the large limit, the radius a should be taken to be significantly 
smaller than 1.9 A, which is the half of the distance between the succesive beads. On the 
other hand, a cannot be too small, since the space along the chain is densely filled with 
amino acids. We take a — 1.5 A in our simulations, which seems a reasonable starting point 
for a qualitative assessment of HI impact on protein unfolding. However, further studies on 
the impact of a on the system dynamics are needed, in particular the hydrodynamic radius 
might need to vary along the chain, reflecting the different sizes of the residues. 

In the approximation Q, the divergence of the diffusion matrix vanishes (Vj • Dy = 0), 
which further simplifies the numerical scheme. However, if the full hydrodynamic interac- 
tions are included, the divergence term should be taken into account!^. 

The simulation using Eq. (J3]) together with ([7]) and ^ will be referred to as Brownian 
Dynamics with hydrodynamic interactions (BDHI) in contrast to a simple BD with the 
diagonal diffusion tensor 

Note that the BD describes configurational evolution of the beads on time scales in which 
the inertia effects of the beads and solvent molecules are negligible^ and, therefore, time 
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scale separation issues discussed in Section 3 are not pertinent here. This feature favours 
BD as a method of choice when simulating stochastically driven motion of proteins at a 
coarse-grained leveP. 

In our previous studies on ubiqutin unfolding * 16 * 42 ^ , we have used LD. Here, on the other 
hand, we incorporate HI within the BD approach. This calls for a comparison of the three 
schemes (LD, BD, and BDHI) to distinguish the effects resulting from HI and those from 
the usage of distinct integration schemes. 

V. CONSTANT VELOCITY STRETCHING 

Fig. 1 presents the force-extension curves for the constant-velocity unfolding at different 
unfolding speeds. In the simulations, both termini of a protein are attached to harmonic 
springs with the elastic constant A;=0.06 e/A 2 . The other end of the N-terminus spring is 
fixed whereas the C-terminus spring moves at a speed v p . We consider three values of v p : 
0.5 , 0.05 and v p = 0.005A/r. 

In the low-speed limit, all the three data sets obtained using the LD, BD and BDHI are 
seen to converge to a single curve. In contrast, at large unfolding speeds, the differences 
between the LD and BD are pronounced. However, strictly speaking, in this time regime BD 
has a limited validity, because of the lack of separation between the momentum relaxation 
time (t v = ™ = 0.5r) and the characteristic time of the aminoacid movement due to the 
stretching t pu u = f- = 3r (for the highest speed quoted above). In the experiments, the 
separation of time scales is huge. In water 7 w 67cr]a k 3 x 10~ 9 g s _1 which leads to 
t v = — « 0.06 ps (for the typical aminoacid mass of m ~ 2 x 10~ 22 g). On the other hand, 
the pulling speeds are of the order of 500 nm/s which gives r p « 0.3 ms, thus there is a 
five-order-of-magnitude separation in time scales. For such a case, LD and BD simulations 
would give exactly the same result. 

The fact that the differences between the BD and BDHI trajectories disappear in the 
limit of small v p is due to the lack of impact of HI at small velocities. To conclude, in the 
experimentally relevant small speed limit, the effects of HI are expected to be negligible. 

An inspection of Fig. 1 indicates that in the case of high stretching speeds neglecting 
HI results in larger peak forces. In the high speed all-atom simulations of titin in water^ 
the forces of stretching are found to be excessively large. Such all-atom molecular dynamics 
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programs are not geared towards hydrodynamic phenomena and may incorporate HI poorly. 
It is possible to consider that the excessive force could be partially due to the missing HI. 

VI. FORCE CLAMP UNFOLDING 

In the force-clamp AFM unfolding j 10 * 11 ! one applies the stretching force to the protein 
terminus and monitors the end-to-end distance, L. The experimental data and the numerical 
simulation ^ 16 * 17 ! show that proteins unfold in a stepwise manner at a constant force. This 
means that a rapid unfolding transition takes place after a certain waiting time. The smaller 
the force, the longer the waiting time. 

In our simulations, we apply the force to the C terminus of the protein, whereas the 
N-terminus is attached to a harmonic spring of elastic constant /c=0.06 e/A 2 . The unfolding 
trajectories of ubiquitin are presented in Fig. 2 for two values of the force, F = 2.4 e/A and 
F = 4 e/A. The LD and BD methods essentially coincide for these relatively large forces 
(small differences between the trajectories are merely stochastic in nature). However, the 
inclusion of HI changes the physics considerably - the waiting times become much smaller 
and the duration of the unfolding transition itself decreases from « 250r to about 50r at 
both values of the force. 

The fact that the HI facilitate protein unfolding may be understood qualitatively when 
one realizes that an amino acid moving away from the bulk of a protein creates a flow which 
drags other residues with it ( see Fig. 3). 

The differences between unfolding with and without HI are further highlighted by analysis 
of the so called unfolding scenarios^, in which one plots an average time when a given contact 
is broken against the contact order, i.e. against the sequential distance, \j — i\, between the 
amino acids that form a native contact. Figs. 4 and 5 compare the unfolding scenarios for 
LD, BD and BDHI at F = 2.4 e/A and F = 4 e/A respectively. Remarkably, although the 
differences in time scales between the unfolding with and without HI are considerable, the 
unfolding scenarios for the smaller force are similar (Fig. 4), which shows that the unfolding 
pathway of a protein is not affected by the hydrodynamic effects. However, as the force is 
increased, both LD and BD scenarios change (Fig. 5): the (3 hairpin structure now unfolds 
at the end instead of at the beginning of the unfolding process. In contrast, in the case 
of BDHI, such a switch is not observed: the scenarios for larger and smaller forces look 
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qualitatively similar. 

VII. UNFOLDING IN A UNIFORM FLOW 

Finally, we study the influence of HI on the characteristics of the protein unfolding in a 
uniform flow with a speed of U. Although the process has not yet been realized experimen- 
tally, the simulation d 42 * 43 * 4 ^ seem to suggest that uniform flow unfolding leads to a richer 
set of metastable conformations than the constant force pulling. For example, when the N 
terminus is anchored, ubiquitin stretches in a flow through two distinct intermediate states 
corresponding to a partial unzipping. 

A detailed analysis of the uniform flow unfolding of ubiqutin in the absence of HI, together 
with the snapshots of intermediate conformations for different anchorings, can be found inpl. 
In that reference, we have mistakenly reported values of the forces in wrong units. Instead 
of the correct unit of e/A we used e/cr, where a was equal to 5A. 

Fig. 6 shows examples of unfolding trajectories of ubiquitin in a uniform flow for four 
different flow velocities, both with and without HI. We observe that unfolding of the system 
with HI requires a much larger flow speed than without. This can be understood qualitatively 
in terms of the so-called no-draining effect^ the residues hidden inside the protein are 
shielded from the flow and thus only a small fraction of the residues experience the full drag 
force of F = — , ~fU (see Fig. 7). In contrast, when no HI are present, this drag force is 
applied to all residues^. 

Notably, although the time scales and velocities involved in the protein unfolding with 
and without HI are completely different, the metastable states are nearly identical (see Fig. 
6). This feature is related to the dynamic character of HI - they do not change the potential 
energy of the system and, therefore, do not affect its stationary properties. In principle, 
however, one could imagine a situation in which, due to the differences in dynamics imposed 
by HI, a system chooses alternative pathways when unfolding. This is clearly not the case 
here which may be related to the directed character of the disturbance (flow in this case or 
the force in force-clamp experiment) which imposes a prefered direction of unfolding, thus 
greatly reducing the set of available unfolding pathways. 

The similarities in unfolding pathways of ubiquitin are further confirmed by the compar- 
ison of unfolding scenarios. As an example, Fig. 8 shows the unfolding scenario for the flow 
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U = 3.5 A/r with HI compared to U — 0.55 A/r without HI (the mean unfolding times are 
comparable in both cases). The scenarios are very close to each other, the main difference 
being that the HI enhance cooperativity by breaking the contacts in a more simultaneous 
fashion. 

VIII. SUMMARY 

In summary, hydrodynamic interactions seem to affect the time scales of unfolding by a 
constant force and by a fluid flow in opposite way but keep the set of the possible metastable 
states. The HI may also reduce peak forces in stretching at a constant speed, although this 
effect weakens with the diminishing stretching speed. 

This work has been supported by the European program IP NaPa through Warsaw Uni- 
versity of Technology. 
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FIGURE CAPTIONS 



Fig. 1. Force-extension curves for the ubiquitin unfolding at a constant speed obtained 
using Langevin Dynamics (dotted line) and Brownian Dynamics with (thick solid 
line) and without (thin solid line) hydrodynamic interactions. The succesive panels 
correspond to the pulling speeds of v p = 0.5,0.05 and 0.005 A/r from top to bottom 
respectively. 

Fig. 2. The end-to-end distance of the model ubiquitin as a function of time during un- 
folding at a constant force. The thick solid, thin solid and dotted lines correspond to 
the BDHI, BD, and LD simulations respectively. The upper panel corresponds to the 
force F = 2.4 e/A, and the lower one to F = 4 e/A. 

Fig. 3. The dragging effect: the moving particle creates a flow patern which affects other 
particles by pulling them in the direction of its motion. 

Fig. 4. Unfolding scenarios of ubiquitin at constant force of F — 2.4 e/A simulated with 
LD (the upper panel), BD (the center panel) and BDHI (the lower panel). Open circles, 
triangles, squares, pentagons and solid triangles and squares correspond to contacts 
(36-44)-(65-72), (12-17)-(23-34), [(l-7),(12-17)]-(65,72), (41-49)-(41-49), (17-27)-(51- 
59), (1-7)— (12-17) respectively. The crosses denote all other contacts. The segment 
(23-34) corresponds to a helix. The two /3-strands (1-7) and ((12-17) form a hairpin. 
The remaining /3-strands are (17-27), (41-49), and (51-59). 

Fig. 5. The same as in Fig. 3 but for F = 4 e/A. 

Fig. 6. The end-to-end distance of the model ubiquitin as a function of time during 
unfolding by a uniform fluid flow as illustrated by several trajectories. The lower 
panel corresponds to the description with the hydrodynamic interactions and the up- 
per - without. The succesive curves in the upper panel correspond to the flows of 
0.25,0.3,0.4 and 0.5 A/r bottom to top respectively. In the lower panel the flow 
speeds are 2.5, 3.5, 5.0, 10.0 A/r In all the simulations, the N-teminus of the protein is 
fixed. 

Fig. 7. The shielding effect: the particles inside a cluster experience a smaller drag force 
than those on the surface. 
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8. The left and right panels show unfolding scenarios without the hydrodynamic 
interactions (U = 0.55 A/r) and with the hydrodynamic interactions (U = 3.5 A/r) 
respectively. 
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